Statistical mechanics of the lattice sphere packing problem 
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We present an efficient Monte Carlo method for the lattice sphere packing problem in d dimensions. 
We use this method to numerically discover de novo the densest lattice sphere packing in dimensions 
9 through 20. Our method goes beyond previous methods not only in exploring higher dimensions 
but also in shedding light on the statistical mechanics underlying the problem in question. We 
observe evidence of a phase transition in the thermodynamic limit d — > oo. In the dimensions 
explored in the present work, the results are consistent with a first-order crystallization transition, 
but leave open the possibility that a glass transition is manifested in higher dimensions. 

PACS numbers: 61.50.Ah, 05.10.Ln, 05.20.Jj 
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The problem of identifying the highest density sphere 
packing in d dimensions is a classical problem in geometry 
with direct connections to fields of physics, information 
theory, and mathematics. The case d = 3, for which 
Kepler conjectured that the face-centered cubic lattice is 
optimal, stood as an open problem for centuries before 
finally in 1998 a proof was announced by Hales Aside 
from dimensions 2 and 3, the highest density is not known 
in any other dimension, although it has been bounded to 
an extremely tight interval in dimensions 8 and 24 [2j. 

The highest packing densities that have been obtained 
in these dimensions, as in many others, are obtained by 
Bravais lattices: periodic packings with one sphere in 
each unit cell. From this point on, we use "lattice" to 
mean a Bravais lattice unless a non-trivial basis is men- 
tioned. When restricting the sphere packing problem to 
lattices, many simplifications are possible, and the prob- 
lem becomes more tractable, but still far from trivial. In 
fact, the densest lattice packings are known in dimensions 
d < 8 and d = 24 0, H| . The space of lattice packings in 
d dimensions is finite-dimensional and much simpler to 
study than the infinite-dimensional space of all possible 
packings. 

In fact, it is possible in principle to identify all lo- 
cal density maxima in this space. Such lattices, called 
extreme lattices, have been characterized by Voronoi in 
terms of their algebraic properties. Voronoi showed that 
a lattice is extreme if and only if it is perfect and eutactic 
Perfect lattices — those that are fully determined by 
a list of their shortest vectors — are finite in number, and 
Voronoi gave an algorithm which enumerates all the per- 
fect lattices in a given dimension. In dimensions d < 8, 
the identity of the densest lattice packing has been es- 
tablished by an exhaustive enumeration of the perfect 
lattices. Voronoi's algorithm relies on a method of ob- 
taining, starting from any perfect lattice, a set of neigh- 
boring perfect lattices. Voronoi showed that this graph 
of neighbors is connected, so exploring larger and larger 
neighborhoods of a single perfect lattice would eventually 
uncover all perfect lattices. 

However, as the number of perfect lattices grows 
rapidly in dimensions d > 8, exhaustive enumeration 
becomes impractical, and other methods must be used 



to identify dense packings. Analytic constructions based 
on groups, codes, and laminated lattices have been suc- 
cessful in producing very dense lattice sphere packings 
in dimensions up to d = 24 and in certain dimensions 
above [3j . While these methods have certainly proved re- 
markably effective, they give little reason, on their own, 
to believe that they have in fact produced the densest 
possible structures. 

For the latter purpose, in the absence of rigorous proofs 
(as in d = 24) or tight bounds, we rely on numeri- 
cal methods that attempt to discover these structures 
de novo: without a priori knowledge of their existence. 
It is only recently that such methods were introduced 
that could tackle moderately high dimensions. A method 
based on the "divide and concur" framework for con- 
straint satisfaction problems was used in Ref. [5| to dis- 
cover de novo the densest known lattices in d < 14. In 
Ref. @, Andreanov and Scardicchio implemented a ran- 
dom walk on Voronoi's graph, yielding a random sample 
of perfect lattices. While their method was designed to 
explore random perfect lattices and their statistics, they 
were also able to use it for de novo discovery of the dens- 
est lattice packing. In dimensions 8 < d < 12, their 
random sample included the densest known lattice pack- 
ing, while in dimensions 13 < d < 19, they had to bias 
their random steps towards higher density configuration 
to recover the densest known packing. In d = 19, only 
some of these biased walks ended up visiting the densest 
known lattice. In Ref. Q, Marcotte and Torquato used 
a sequential linear programming (SLP) approach to iter- 
atively compress a lattice configuration until it reaches a 
configuration which cannot be compressed any more — an 
extreme lattice. This procedure, starting from random 
initial conditions, reproduced the densest known lattice 
in an appreciable portion of the runs in each dimension 
d < 16. The percentage of runs yielding the densest 
known lattice declines sharply starting at d — 17. How- 
ever, because each run can still be computed rapidly, the 
procedure can be repeated many times and the densest 
known lattice is produced at a decent rate for d < 19 (see 
below for a direct comparison with the present method). 
Other methods have been used for de novo searches in 
closely related problems, such as the Gaussian-core soft 
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sphere ground state problem [8| and the lattice quantizer 
problem [jj. 

In this paper we report on a Monte Carlo (MC) method 
for studying the lattice sphere packing problem. Our 
method is completely different from the references above, 
but some of our results are surprisingly similar to the re- 
sults of Refs. 0, 0]- In particular, in each dimension 
d < 19, a simulated quasistatic comression discovers de 
novo the densest known lattice in at least 30% of the runs. 
In 20 dimensions only one out of 50 quasistatic compres- 
sion runs yielded the densest known lattice. However, 
with a slight change in protocol, we were able to repro- 
duce the lattice in 7 out of 50 runs. That three such 
disparate methods all seem to start having serious diffi- 
culties in the same dimensions raises the possibility that 
the lattice sphere packing problem becomes intrinsically 
harder around d — 20. We suggest possible reasons for 
such a scenario. Moreover, because of the statistical me- 
chanical nature of our method, we are able to quanti- 
tatively characterize the intrinsic nature of the lattice 
sphere packing problem apart from the behavior of any 
specific method or algorithm for its solution. 

A lattice can be specified by its generating matrix: the 
set of all sphere centers is given by M Tr L d = {M T n : n £ 
Z d }, where M is a d X d matrix and its rows are gen- 
erating vectors (primitive vectors) of the lattice. While 
a generating matrix uniquely specifies a lattice, a lattice 
has multiple generating matrices: whenever Q is a uni- 
modular integer matrix, M' = QM generates the same 
lattice as M. A lattice is a packing of radius- 1/2 spheres 
(that is, its spheres do not overlap) if ||M T n|| 2 > 1 for 
all n G 7L d \ {0}. If this is the case, we say the lattice 
is admissible. The number density of spheres is given by 
1/v = l/|detM|, where v denotes the unit cell volume. 
Any lattice, generated by M, can be rotated so that its 
generating matrix MU, where U is a rotation matrix, 
becomes lower-triangular. Therefore, the space of lat- 
tices, modulo rotation, can be parameterized by lower- 
triangular generating matrices. 

We can define an isobaric ensemble on the space of 
admissible lattices by weighting the probability of each 
lattice by a factor exp(— pv), where p = NP/ksT is a 
reduced pressure variable (cf. Ref. 10] for definitions 
of related ensembles). If v is thought of as the energy 
associated with a particular lattice, we can think of p as 
the inverse-temperature (cf. Ref. [6|, where log?; is used 
as the energy). We can sample this ensemble using a 
standard Metropolis algorithm: a random element of the 
lower-triangular generating matrix is randomly changed 
by a small amount; the step is rejected always if it yields 
an inadmissible lattice and is rejected with probability 
exp(— pAv) if the unit cell volume is increased by Av. 
However, we find instead that it is more efficient to build 
the detailed balance directly into the proposed steps in- 
stead of the acceptance probability. Note that changes 
to the off-diagonal elements of the matrix do not change 
the volume. Therefore, these moves are always accepted 
if they produce admissible lattices. For the diagonal ele- 



ments, we propose changes 

/ ex-e 2 p\ 
TTla 4 IH 1 ma, 

where e is a measure of the typical move size and x 
drawn from a normal distribution. As with the off- 
diagonal moves, the proposed moves are always accepted 
if they produce admissible lattices. Note that the prob- 
ability of a change Av in the volume is proportional to 
exp((Au — e 2 p)/2e 2 ), so an accepted move changing the 
volume by Av is more likely by a factor of exp(pAv) 
than the reverse move. This bias is the only place that 
the pressure enters in our method. 

A crucial step in this MC algorithm is checking 
whether a lattice is admissible. This is known to be an 
NP-complete problem [llj, and in fact it takes up most 
of the computational time in our simulations. The com- 
plexity of this problem is sensitive to the choice of gener- 
ating matrix for a given lattice. Generally speaking, the 
shorter and more orthogonal to each other the generating 
vectors are, the easier the problem is of determining ad- 
missibility of the lattice they generate. There are many 
non-equivalent criteria for determining how well-suited, 
or reduced, a certain set of generating vectors is. We use 
Korkine-Zolotareff (KZ) reduction, which is one of the 
most stringent of these criteria [l2[ • Such a stringent cri- 
terion is warranted because of the large number of times 
we are required to decide the admissibility of similar lat- 
tices. Therefore, during our simulation we periodically 
perform KZ reduction, and thus all the generating ma- 
trices we consider are either KZ-reduced or nearly so. 
Detailed balanced is violated by these reductions, but we 
find that the reductions are rare enough in the key stages 
of the simulations that they do not significantly hinder 
thermalization. 

Using our MC technique, we perform a simulated qua- 
sistatic compression (a simulated annealing where pres- 
sure takes the role of temperature) . We start the system 
in a simple hypercubic lattice Z d and equilibrate at a 
constant pressure. We then start to increase the pres- 
sure by a constant factor after each proposed move. We 
vary the typical move size inversely with the pressure: 
e = e /p. We use different move sizes for off-diagonal 
and diagonal moves, and for both we pick eo so as to 
achieve an average move acceptance rate of roughly 30%. 
The length of the equilibration period is 4% of the length 
of the compression period. The Gram matrix G — MM T 
of the final matrix obtained in nearly all runs consists, up 
to small errors, of small-denominator rational numbers. 
Therefore, we can easily round off its entries to obtain 
the infinite pressure limit of the simulation. 

For each dimension d = 9, . . . , 19 we performed 20 
simulation runs. In each of these dimension, our sim- 
ulations discover the densest known lattice packing in at 
least 30% of the runs. Table [J summarizes the results 
of these simulations. The compression rate used repre- 
sents a trade-off between longer computation time and 
decreased likelihood of reproducing the densest lattice. 
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FIG. 1: Traces of density as a function of reduced pressure 
in different runs in dimensions d = 11 and d — 16. In 11 di- 
mensions, the simulation remains in equilibrium until around 
p = 2.4 x 10 3 , where different runs continue along different 
branches, corresponding to basins of attractions of different 
extreme lattices. In 16 dimensions, all runs end up in the 
basin of the same lattice, but we observe that different runs 
make the transition into this basin at different pressures. It 
appears that at least one of the runs also spends some time 
in an intermediate state, presumably the basin of attraction 
of a different extreme lattice. 



For any fixed d, as there is no thermodynamic limit, 
strictly speaking, there cannot be a phase transition. 
However, as is clear from Figures [T] and [2j the system 
shifts as the pressure increases from a state where many 
basins of attraction are explored to a state where the sys- 
tem is confined to a single basin. In any finite dimension 
there should be a range of pressure where these two states 
coexist with significant probability for the system to be 
in either state. The traces we obtain from the simula- 
tions are consistent with a situation where the transition 
rate between these two states in the coexistence region 
becomes smaller and smaller with increased dimension, 
so that in lower dimensions the trace of each run remains 
close to the equation of state, while in higher dimensions 
each run stays in one state until transitioning irreversibly 
into the other. This crystallization transition is accom- 
panied by a discontinuity in the density, which depends 
on the extreme lattice the system crystallizes into. 

We may interpret d — > 00 as the thermodynamic limit 
of the lattice sphere packing problem and speculate that 
the coexistence region in this limit shrinks to a single crit- 



TABLE I: For each dimension d, the table gives the follow- 
ing: the name (as per Refs. % Q]| ) of the lattice A that 
achieves the greatest density known among all admissible lat- 
tices; the unit cell volume v of this lattice (normalized by 2 d ); 
the reduced pressure pi used in the equilibration period and 
at the beginning of the compression period of our simulations; 
the reduced pressure Pf at which we terminated the compres- 
sion; the compression rate k, such that the pressure at each 
proposed move is 1 + k times the pressure at the previous pro- 
posed move; the average computational time T per run of the 
simulation; and the rate at which the lattice A is reproduced, 
that is, the percentage of runs whose final configuration, after 
rounding off the Gram matrix, is A. 



We did not attempt to quantify this trade-off in this pa- 
per or determine the optimal compression rate as a func- 
tion of dimension. In terms of the average computational 
time needed to reproduce the densest known lattice once, 
the present method is much less efficient than the SLP 
method of Ref. Q in all lower dimensions. For the high- 
est dimensions, d = 17, 18, and 19 resp., this average time 
is 3 x 10 3 s, 8 x 10 4 s, and 4 x 10 6 s with the SLP method, 
compared to 2 x 10 5 s, 3 x 10 5 s, and 8 x 10 5 s with the 
present method. 

Assuming our MC simulations accurately sampled the 
isobaric ensemble at each pressure, we would recover from 
the simulation the equation of state: the average volume 
(v) as a function of the reduced pressure p. Plotting the 
traces of (v) as a function of p for different runs in di- 
mension d = 11 for example (Figure [T]), we see that at a 
certain pressure the different traces diverge and the sim- 
ulations fall out of equilibrium. The traces belonging to 
runs that terminate at the same configuration do not di- 
verge, so we may infer that at this pressure the system 
goes from exploring the attraction basins of many dif- 
ferent extreme lattices to exploring only a single basin. 
The situation becomes more complicated in higher di- 
mensions, where we see significant hysteresis effects. For 
example, in dimension d = 16, all 20 runs end up in the 
basin of the densest known lattice Aie (Figure [I}, but 
different runs experience the transition at different pres- 
sures. Figure [5] shows, for each dimension, the volume 
as a function of pressure averaged over all the runs that 
yielded the densest known packing. 
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FIG. 2: Traces of (normalized) density as a function of re- 
duced pressure, averaged over all quasistatic compression runs 
that yielded the densest known lattice in a each dimension 
d = 9, ...,20. 



ical pressure. Parisi considers the partition function of a 
closely related ensemble of lattices in the limit d oo 
and points out that the lattice sphere packing problem in 
this limit shares many formal similarities with the non- 
lattice problem [l(| . Parisi does not determine whether a 
glass transition exists in the lattice sphere packing prob- 
lem as it does in non-lattice hard sphere system, and 
leaves open the possibility of either a glass transition or 
a crystallization transition. The behavior we observe in 
our simulations in dimensions 9 < d < 19 is indicative 
of a first-order crystallization transition. However, it is 
hard to tell whether the behavior is controlled by the 
thermodynamic limit or mostly by details specific to the 
moderate dimensions we explore. 

In dimension d — 20, out of 50 runs at the slowest 
compression rate attempted, only one yielded the densest 
known lattice. In fact, only ten runs show a discontinuity 
in the density at all, with most runs remaining in the 
fluid state throughout. In another set of 50 runs, we 
compressed to an intermediate pressure, where we expect 
the crystallization rate to be higher, and maintained that 
pressure until we observed a rapid increase in density. 
With this new protocol, the densest lattice is reproduced 
in seven runs. 

We find it remarkable that both our method and the 
methods of Refs. d, 0| become dramatically less effective 
at exactly the same dimension. Whereas a priori, we 
might expect that the complexity of the lattice sphere 
packing problem rises at a more or less constant expo- 
nential rate as a function of dimension, the evidence of 
the two cited works and the present paper raises the hy- 
pothesis that the complexity experiences a sharp increase 
around d = 20. A sharp increase of this kind might in- 
dicate a shift into the glassy regime of the lattice sphere 
packing problem. Just as the relaxation time of a fluid 
increases sharply as the glass temperature is crossed, it 
might be the case here that d « 20 marks the begin- 
ning of a glassy regime, linked to a sharp increase in the 
inverse compression rate required to recover the densest 



lattice. 

In addition to discovering de novo the densest known 
lattice in dimension d, our method can also be used to 
discover suboptimal, yet very dense, extreme lattices. In 
some dimensions, only a portion of runs, even at the 
slowest compression rate attempted, yielded the densest 
known lattices (see Table H}. The identity of the sub- 
optimal lattices produced is in some cases unexpected, 
and the second-most-likely-produced lattice is not always 
the second-densest known extreme lattice. For exam- 
ple, in dimension d = 14, the most frequently produced 
lattice after A14 in our simulations is a lattice (denoted 
"diml4kis744" in Ref. [l3[) of normalized unit cell vol- 
ume 2 14 v = 361-^/3/16, compared to 2 14 v — 16\/3 for 
A 14 . This is also the second-densest lattice produced, 
despite the existence of many extreme lattice of interme- 
diate density [7|. As was already observed in the results 
of Ref. Q , we observe a general trend wherein among ex- 
treme lattices of equal density, those with lower kissing 
numbers (number of neighbors in the first coordination 
shell) are more frequently produced, though this trend is 
not without exceptions. In the present context, it makes 
sense that at finite pressures the basin of an extreme lat- 
tice with a lower kissing number is stabilized by a greater 
rattling entropy over the basin of another extreme lattice 
of equal density and higher kissing number. 

Many of the lattices discovered by our simulations have 
not, to our knowledge, been studied before, and we have 
submitted them to be archived in the on-line catalog of 
lattices [l3j]. A few in particular are definitely worthy 
of further study. For example, one of the lattices we 
discover in dimension d — 11 (denoted "dimllkis422" in 
Ref. [l3|]) is equal in density to the two laminated lattices 
A™ m ' max . Remarkably, the lattice does not include any 
of the laminated lattices A^ for d — 8, 9, 10 as sublattices 
of equal minimum norm. This discovery suggests a pos- 
sible extension to the conventional lamination hierarchy 
described in Refs. HEIEU 

While the focus of the present paper is limited to the 
lattice sphere packing problem, the method presented can 
easily and naturally be extended to study lattices with 
an n-element basis for n > 1. The smallest dimension in 
which the densest known packing has a non-trivial ba- 
sis is d = 10, where it has a 40-element basis Q. In 
higher dimensions, there are known packings with more 
modest basis sizes (for example n = 4, 3 in dimensions 
d = 20,22 resp. [ill) that are denser than the densest 
known lattices. Extending our capabilities to the point 
of being able to discover de novo any of these non-lattice 
structures or the densest known lattices up to d = 24, 
would be a major accomplishment. Of course, discover- 
ing yet unknown lattice and non-lattice packings denser 
than those constructed analytically should be considered 
the ultimate goal. 
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